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(N ; ABSTRACT 

Qh 1 In the second paper of this series we pursue two objectives. First, in order 

to make the code more sensitive to small effects, we remove many approxima- 
■ tions made in Paper I. Second, we include turbulence and rotation in the two- 

dimensional framework. The stellar equilibrium is described by means of a set 
of five differential equations, with the introduction of a new dependent variable, 
namely the perturbation to the radial gravity, that is found when the non-radial 
effects are considered in the solution of the Poisson equation; following the scheme 
of the first paper, we write the equations in such a way that the two-dimensional 
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effects can be easily disentangled. The key concept introduced in this series 
is the equipotential surface. We use the underlying cause-effect relation to de- 



velop a recurrence relation to calculate the equipotential surface functions for 
uniform rotation, differential rotation, rotation-like toroidal magnetic fields and 



turbulence. We also develop a more precise code to numerically solve the two- 
dimensional stellar structure and evolution equations based on the equipotential 
surface calculations. We have shown that with this formulation we can achieve 
the precision required by observations by appropriately selecting the convergence 
criterion. Several examples are presented to show that the method works well. 
Since we are interested in modeling the effects of a dynamo-type field on the 
detailed envelope structure and global properties of the Sun, the code has been 
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optimized for short timescales phenomena (down to 1 yr). The time dependence 
of the code has so far been tested exclusively to address such problems. 

Subject headings: Sun: evolution — Sun: interior — stars: variables: other - 
Sun: Oscillations 



1. INTRODUCTION 

High precision is an essential requirement in solar variability modeling because the 
cyclical variations of all solar global parameters are very small (see Li et al 2003 and references 
therein). For example, the (relative) precision of the measurements of the TSI is about 10~ 5 . 
Oscillation splittings can also be measured with a similar precision, and the PICARD satellite 
expects to measure diameter changes with a precision of a few milli-arc seconds, thus a few 
parts in 10 6 . These requirements are even more extreme in the two-dimensional (2D) case, 
because two-dimensional effects are subtler than their ID counterparts. This gives us a sense 
of the precision required for our code. 

In the first paper of this series (Li et al 2006, referred hereafter as Paper I), we developed 
a 2D stellar evolution code that includes magnetic fields of arbitrary cylindrically symmetric 
configuration by generalizing in a straightforward way our one- dimensional (ID) code (Lydon 
& Sofia 1995; Li & Sofia 2001; Li et al 2002; Li et al 2003). Since the 2D case is very complex, 
we made some significant approximations, some physical, and some computational. In terms 
of the physical approximations, for the first two, we neglected the second-order derivative 
of the gravitational potential $ with respect to the colatitude coordinate 9 and the second- 
order derivative of the perturbation gravitational potential $ — $ with respect to the radial 
coordinate r, where d^o/dr = Gm/r 2 is the spherically-symmetric gravitation acceleration 
component in the radial direction, i.e., expression (30) in Paper I. The third approximation 
is that we ignored turbulence, which had been included in our ID variability models (Li 
et al 2002). A detailed comparison of the ID solar variability models with the relevant 
observations (Li et al 2003) shows that turbulence must play an important role. In particular, 
in order to explain the changes of the oscillation spectrum in function of the activity cycle, 
we needed to include a model of turbulence that interacts with magnetic fields in a negative 
feedback sense. In this paper we remove these three physical approximations made in Paper 
I. 

Unlike the three approximations mentioned above, the fourth approximation made in 
Paper I is computational, involving the solution method of the 2D stellar structure equations. 
In the ID case, we use the trapezoidal rule to integrate the ID stellar structure equations. 
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In the 2D version in Paper I, the trapezoidal rule (or the central difference scheme) was 
not applied everywhere, since we used numerical derivatives. In this paper we minimize the 
use of numerical derivatives. The fifth approximation made in Paper I is that we neglected 
OFq/BO in the luminosity equation, i.e., the term 0(2) in Eq. (124d), which we now include. 
The similar term in Eq. (124e) of Paper I does not matter for the cyclic variation of the Sun. 

Removal of the above six approximations is one of the main objectives of this paper. 
The second main objective is to include turbulence and rotation, which are also important 
sources for asphericity. In §2] we summarize the theoretical foundations that give rise to the 
2D stellar variability models by including magnetic fields, turbulence and rotation. Since 
we want to get rid of approximations 1 and 2, we have to add the Poisson equation (which 
is a second-order partial differential equation) to the stellar structure equations. We thus 
have two more first-order stellar structure equations. As a result, we now have a total of six 
stellar structure equations. 

Equipotential surface is the key concept to obtain the 2D generalization from the ID 
stellar structure and evolution equations. In §3] we show how to find out the equipotential 
surface from the 2D stellar structure equations obtained in §2J Magnetic fields, turbulence 
and rotation are causes, and the resultant matter redistribution is the effect. This cause- 
effect relation indicates certain recurrence relation for equipotential surface calculations. We 
present the recurrence relations for the uniform rotation, differential rotations, rotation-like 
toroidal magnetic fields and turbulence in §3J 

The third main objective of this paper is to raise the numerical precision of the numerical 
solutions for the 2D stellar structure and evolution equations. We tried hard to do so and 
found out that it is the best to explicitly invoke the equipotential surface. We present this 
method of solution in §H We give a typical example of 2D solar variability models in §5] to 
show how we use this 2D code. The conclusion is presented in the last section. 

2. THEORETICAL FOUNDATION 

Magnetic fields, turbulence, and rotation are possible causes of asphericity. In this 
paper, we consider all of them. We assume that the system is azimuthally symmetric or 
axisymmetric. Therefore, we need only the radius (r) and colatitude (9) in the spherical 
polar coordinate (r,9,(p), in which the azimuthal angle <fi is irrelevant. The basic equations 
represent the conservation of mass, momentum, and energy. We also need the Poisson 
equation and the energy transport equation to close the system. Since magnetic fields are 
involved, the Maxwell equations must also be obeyed, for example, we require V ■ B = 0. In 
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this section we summarize the results and point out the differences from their ID counterpart. 



2.1. Mass Conservation 

Mass conservation is guaranteed by calculating the mass enclosed within a certain sur- 
face. In a spherically symmetric system the surface is a spherical surface with radius r with 
respect to the symmetric center of the system. This spherical surface is also an equipotential 
surface of gravity. In the general case the equipotential surface r = 9) is thus used to 
define the mass M$: 

M$ = 2n d9sm9 drr 2 p(r,9). (1) 

Jo Jo 

In the spherically symmetric case we have M$ = M r since the equipotential surface R(<&,9) 
and the density p(r, 9) do not depend upon colatitude 9. 

This mass expression sets up a one-to-one relationship between mass M$ and equipo- 
tential $: 

m = M$ = M$($), (2) 

which permits us to use mass M$ and colatitude 9 as the 2D independent variables, instead 
of the gravitational potential $ and colatitude 9. Eq. (TjQ) is the integral form of the mass 
conservation. Its differential form can be obtained by taking its partial derivative with 
respect to the equipotential surface r = 9): 

dm dMs, , 9 

= — — ATrr z 

dr dR 

where 



47rr p m , (3) 



Pm = — / d9R 2 ($,9)p(R($,9),9)sm9. (4) 
Zr Jo 

This defines the density on the equipotential surface r = i?($,6 l ). It should be pointed out 
that here r is no longer a static Eulerian space coordinate, but a co-moving Lagrangian 
variable with an equipotential surface $ = constant, as it is in the spherically symmetric 
case. Obviously, p m = p when the system is spherically symmetric. We use r = r(m,9) to 
denote the functional relationship between r and m. 

Eq. d3J) is our 2D mass conservation equation. Comparing it with its ID counterpart, 

dm 9 . , 

— = Anr 2 p, (5) 



we find that they differ by a correction factor p m j p: 

dm 2 j p, 

_ = 47IT pi — 

or {p 



dm -4.r 2 pl^}. (6) 
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This factor equals unity in the spherically symmetric case, but deviates from unity in the 
general case. 



2.2. Momentum Conservation 

When both turbulence and magnetic fields are taken into account, the momentum con- 
servation of an equilibrium state can be expressed by the momentum equation 



P + f- + P« ) 1 + pWe ~ «)eee e 



in 



+PW - v'X)e^ - i-BB 



-pV$ - V • (pvv), (7) 



where P is the gas pressure, B is the magnetic field, I is the unit tensor with nonzero 
components e r e r , e d e d and e^e^,, and v" is the turbulent velocity that is defined by the 
velocity variance: 



< = ^-vff'\ (8) 

where i — r,9, 0. The over-bar denotes a combined horizontal and temporal average, and t>j 
is the total velocity component. See Robinson et al (2003) for the details of 3D simulations to 
derive realistic turbulence properties in the solar convection zone, where v' e ' = v'i is assumed. 
The regular motion velocity is denoted by v, for example, v = Q x r for rotation, where ft 
is the rotation angular velocity. 

For a system with magnetic fields, turbulence and rotation, Eq. (jTJ) can be rewritten as 
follows: 

-jr- = -p^- + p(n r + % + n r ), (9) 

or or 
ldP r pd§ 

~ = --M+P(H + T e + Il e ), (10) 
r do r do 

where the isotropic pressure components of the magnetic field B, P m = B 2 /8n, and the 
radial pressure component of turbulence, P t = pv"v" have been added to the gas pressure, 
P, to define a total isotropic pressure Pt = P + P m + Pt, while their anisotropic pressure 
components are denoted by TC = V ■ (BB) for the magnetic field B, T = p _1 V ■ [p(v"v" — 
v'gv'^egeg + p{v"v" — v'Lv'^e^e^ for turbulence, and 71 = — V ■ (pvv) for rotation, where 
v = Q x r. Their r- and ^-components are: 

^ P H r = ^yBl) + \^{B r B 6 )-\{Bl + Bl), (11a) 
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n 2 r sin 2 6, 
VL 2 r sin 6 1 cos ^. 



cot 9, 



(lib) 
(11c) 

(lid) 

(lie) 
(llf) 



2.3. Poisson Equation 

The Poisson equation in the spherical coordinate system with the specified symmetry 
requirement can be written down as follows : 

1 d ( 2 <9$\ 1 d ( . „d<$> 



r 2 dr \ dr ) ~^ r 2 sin 9 89 \ 39 ) ^ ^ ^ 

Solving this equation for the gravitational potential $ is not sufficiently accurate for our 
purposes, especially in the core of stars. Solving it for the radial gravitational acceleration 
g = d^/dr is equally not good for the same reason. Many tries show that the following 
treatment is sufficiently accurate for our high-precision requirement. 

First of all, we calculate the colatitudinal gravitational acceleration Q = (l/r)d&/d9 by 
using the hydrostatic equilibrium equation in the colatitudinal direction (Eq. [TO]) in terms 
of dP T /de, H e , and T e : 

1 dP T 

g = H e + Te + n e -f. (13) 

rp o9 

This way, Eq. ffTUl) is satisfied automatically. We then decompose g into two parts, 

+ Sg. (14) 



,2 



The first part is the spherically-symmetric radial component of the gravitational accelera- 
tion, and the second part is the deviation of the radial gravitational acceleration from its 
spherically-symmetric counterpart. Substituting Eq. (fT4l) into Eq. (fT2"l) . we obtain 



(15) 
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Therefore, we solve the Poisson equation for Sg instead of $ or g. Here we have used 
the notations r' = lnr, and p' = In p. The hydrostatic equilibrium equation in the radial 
direction thus becomes: 

^ = -p(^ + 5g-H r -T r -n\ (16) 
or \ r z I 



2.4. Energy Conservation 

The energy conservation equation is 



r 2 dr^ r ^ r sin 6 39^' P y dt j ^ ^ 

where F = F ra d + F conv is the energy flux vector, including both the radiative flux F ra d and 
the convective flux F conv , and e is the rate of nuclear energy generation, and St is the total 
specific entropy, including the contributions from magnetic fields and turbulence. We use 
the diffusion approximation for radiative flux, and the mixing length theory for convective 
flux: 

AacT 3 



rad 



-VT, 



F conv = _1 Pj^conv (lg) 

21 + v conv / v Q 

where f con v is the convection velocity, l m is the mixing length, v is a typical velocity deter- 
mined by choice of radiative loss mechanism of a convective eddy. The symbol a represents 
the radiation constant, c the speed of light, k the mass opacity coefficient. The 2D en- 
ergy conservation equation shows that energy can not only penetrate a region via the radial 
gradient of the radial component of the energy flux, but also goes around it via the trans- 
verse gradient of the transverse component of the energy flux. In contrast, the ID energy 
conservation equation 

1 d , 2 . ( dS T 



--y Fr) = p U-T^) (20) 



rules out the transverse transport of energy. 



2.5. Energy Transport 



Eqs. (I17ti20p show that we have to calculate temperature and entropy gradients. We 
thus need the first law of thermodynamics in the presence of magnetic fields and turbulence. 
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We have redefined the mechanical variable Pt by adding all isotropic pressure components 
together. We need magnetic and turbulent variables to take into account magnetic and 
turbulent degrees of freedom. 



2.5.1. Magnetic and turbulent variables 

We use B to define three stellar magnetic parameters, in addition to the conventional 
stellar parameters such as pressure, temperature, radius and luminosity. The first magnetic 
parameter is the magnetic kinetic energy per unit mass, Xm, 

Xm = B 2 /(87ip). (21) 

The second is the heat index due to the magnetic field, or the ratio of the magnetic pressure 
in the radial direction to the magnetic energy density, 7 m — 1, 

lm = l + {B 2 e +Bl)/B\ (22) 

The third one is the ratio of the magnetic pressure in the colatitude direction to the magnetic 
energy density, -# m — 1, 

$ m = l + {Bl + Bl)/B 2 . (23) 

We can use these three magnetic parameters to express three components of a magnetic field 
as follows: 

B r = [87r(2- 7m ) Xm p] 1/2 , (24a) 
B e = [87r(2-^ m ) Xm p] 1 / 2 , (24b) 
B4, = [8vr( 7m + ^ m -3)x m p] 1/2 . (24c) 

However, since v'g = v'^ is assumed, we have only two turbulent degrees of freedom and we 
thus need two turbulent variables, namely, the turbulent kinetic energy per unit mass, Xt, 
and the effective ratio of specific heats due to turbulence, 7 t : 

Xt = \{v")\ lt = l + 2{v';/v"f- (25) 
We can use them to express three turbulent velocity components: 



< = [(lt-l)Xt] 1/2 , (26a) 

n l/2 

(26b) 



v'l = v'i 
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2.5.2. Equation of state 

Using the magnetic and turbulent variables denned above, we can rewrite the total 
pressure as follows: 

P T = P(p, T) + p Xm + p{lt - l)Xt- (27) 

Solving this equation for p, we obtain the equation of state in the presence of magnetic fields 
and turbulence: 

p = p(Pr,T,Xm,Xt,lt)- (28) 

To highlight magnetic and turbulence effects we adopt a given chemical composition. This 
shows that the independent thermodynamical variables are Pt, T, Xm, Xt, and j t - Using 
them,we can write the differential form of the equation of state as follows: 

dp/p = adP T /P T - 5dT/T - v m dXm/Xm ~ v t dxt/xt ~ fitdjt/jt, (29) 

where 

a = (dlnp/d\nP T ) T>XmiXt)lt , 5 = -(<91np/<91nT)p TiXmiXti7t , (30a) 

u m = -(dlnp/dlnXm)p T ,T, Xt , lt , u t = -(d\n p/d\nxt)p T ,T, Xm ,~, t (30b) 

fit = -(<91np/<9hi7 t )p TiTiXmiXt . (30c) 

When a ^-dependent magnetic field is applied, Eq. (T28"|) demonstrates that the mass 
distribution will adjust to generate asphericity. This is the most straightforward 2D effect. 

2.5.3. The first law of thermodynamics in the presence of magnetic fields and turbulence 

The first law of thermodynamics is the energy transfer and conservation law in a ther- 
modynamic system. In the presence of magnetic fields and turbulence, the conservation law 
should be modified as follows: 

TdS T = dU + PdV -dxm-dxt, (31) 

which states that both magnetic and turbulent energy are generated at the expense of internal 
energy of the system U. Here V = 1/p is the specific volume. Combining Eqs. (|28|) and (13TJ) 
(see Lydon & Sofia 1995 for the detail), we obtain 

TdS T = C p dT -( 6 -)dP T +( ^ -l)d Xm +(^-l) d Xt + (32) 
\Pj \ctpXm J \apxt J ctplt 
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from which we obtain 



cISt 



{C p /T)VT-{C p VJP t )VPt, 

(c P /T) d S- { c P vjp T ) d -^. 



We have defined the modified adiabatic gradient 



ad 



ad 



1 - I — 

a 



C P T 



a 



Xt 

C P T 



a 



(33) 
(34) 

(35) 



where C p is the specific heat per unit mass at constant total pressure, constant magnetic 
energy per unit mass, constant turbulent kinetic energy per unit mass, and constant turbulent 
specific heat ratio, and 



V ad = P T 5/( P C p T) } V r 



<91nP T ' dlnPr Wy d\nP T 

The physical meaning of Eq. (1331) is that magnetic fields and turbulence provide additional 
channels for energy transport. 



2.5.4- Energy flux vector 

Using Eqs. (I18til9j) and fl33l) the energy flux vector F can be expressed by the temperature 
gradient VT and pressure gradient VPt as follows: 

F _ / 4acT 3 1 pC P l m v conv \ ^ T 
V 3k P 2 1 + v COQV /v J 
l pCpTV'J m v com J_ 
+ 2 l + v conY /v Pt T ' K } 

Its r-component determines the radial temperature gradient, its ^-component results in 2D 
effect. 



3. EQUIPOTENTIAL SURFACE 

Solar magnetic fields are weak in the sense that the resultant magnetic pressure is much 
smaller than the gas pressure. The usual central difference scheme alone may not discern 
the required 2D effects. We should therefore use certain physical guidelines to improve the 
precision of the numerical solutions for the 2D stellar structure equations. The key concept 
introduced for the 2D stellar structure in this series is the equipotential surface. In this 
section we show how to determine it. 
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3.1. Exact 2D Stellar Structure Equations 



The exact 2D stellar structure equations, i.e., Eqs. (J6J), (1131) . ffT5]) . f)T6l) . ffTTj) . and the 
energy transport equation, can be rewritten as follows after coordinate transformation from 
(r,6) to (m,6): 



dr' 
ds 

dp' 

ds 

dr 

ds 

dL 

ds 

dU 

ds 



m p 
4:7ir 3 p p m ' 

m p ( Gm rT ^_ ^ 

dP' [ V ra d radiative 
V c convective 

, dSr \ p 1 ttiFq cot 9 

dt J p m 

m 



ds 
i 

Gm ( p 



-m e 



T- 



1 m dFp 



- 1 



Pr, 



Aur^pr, 



Lq rp m 

2U + g cot e + 



L Q rp m d6 

eg; 
de 



(37a) 
(37b) 

(37c) 

(37d) 
(37e) 



Here P' = lnP^, T' = InT, r' — lnr, L — 47rr 2 F r / L Q , and U = 5g. The other symbols used 
above are defined as follows: 



Q 



AacT 4 + 1 pC P Tl m v, 
3np 

% + TZf, 



V , lpC P Tl m v conv V' ad \ dP' 



2 1+ v com /v J r 2 1 + v com /v a r 
P T dP' 
rp d9 



de ' 



(38a) 
(38b) 



These equations show that in addition to the dependent variables, pressure Pt, temper- 
ature T, radius r, and luminosity L, we have two more dependent variables, the radial and 
colatitudinal gravitational acceleration perturbations 5g and Q. However, we need to solve 
only five partial differential equations (Eqs. l37aH37el) because the colatitudinal gravitational 
acceleration Q can be calculated by using , {%■) > 'He, %, and TZq. 

We use Sg = at m = as the central boundary condition for the fifth equation because 
Sg is a perturbation in nature. This is equivalent to assume that the radial gravitational 
acceleration be equal to its spherically symmetric counterpart at the center. 



3.2. Equipotential Surface Profile 



We know that r is the radial coordinate of an equipotential surface. Its dependence on 
the colatitudinal coordinate 9, i.e., r = r($, 9) = r(m,9) defines an equipotential surface 
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on which the potential equals $. We redefine r by r e = r e (m) and x = x(m,9): r = 
r e (m)x(m, 9), where r e is the equatorial radius. Since we interpret r e as the equatorial radius, 
x should always be normalized so that we obtain x — 1 at the equator where 9 = 7r/2. 
The equipotential surface is thus expressed by x = x(m,9), which is a function of mass 
m = M$($) and colatitude 9. 

In order to find out the equipotential surface x, we use the fact that pressure is 9- 
independent on it. Otherwise, the hydrostatic equilibrium is not reached thereon. This 
indicates that should be ^-independent thereon as well. The following equation is 9- 
independent and holds well for both spherically-symmetric and aspherical cases: 

dP' Cm 2 

ds Aixr A e P T v ; 

Comparing it with Eq. f!37bj) . we obtain 

A P 



x 

Pr, 



r 2 2 



1 + -77- ( u - n r - % - Tl r 

Gm 



(40) 



We can use the iteration method to obtain x starting from x = 1 if we know r e , U, 7i r , %■ 
and 7Z r . Eq. (jlj) that is used to determine p/p m now becomes 

I /-t/2 



Pm = I px 2 sin 9d9. (41) 
x 2 



3.2.1. Mass conservation for r e 

To calculate r e , we rewrite the mass conservation equation as follows: 

dr 

^ = 1/Q ' (42 > 

where 

Q = ATtr 2 Pm (43) 

is independent. As a result, we know that J^- is ^-independent. Therefore, we can choose 
r at any specific colatitude on the equipotential surface. We can, of course, choose r = r e (m) 
to obtain 

dr>„ m 



e 



ds Qr e 

where ri = lnr P . 



e 



Eq. (1441) becomes 

m = -Qr e (45) 
at the center. This is one of the central boundary conditions. 
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3.2.2. Poisson equation for U 

The radial gravitational acceleration perturbation U = 5g can be decomposed into five 
components U = Ud + Up + Uh + Ut + Ur according to their physical origins specified by the 
subscripts, where subscript D stands for the density variation, P for the pressure variation, 
H for magnetic fields, T for turbulence, and R for rotation. To see this, we decompose the 
colatitudinal gravitational acceleration component into four components according to their 
physical ingredients G = Gp + Gh + Gt + Gr- Their definitions are 

GmQ ( dx' 



Anrip V M ' (46a) 

Gh = He, (46b) 

Gt = %, (46c) 

Gr = Tie, (46d) 

where we have utilized the equipotential surface condition (^-) = and defined x' = In a;. 
Since the Poisson equation is linear, we can write it down for each component as follows: 

-rr 1 = + Si, (47) 

or r 

where i = D, P, H, T, and R. The source terms Si are expressed by the following functions: 

S D = 4irG(p- Pm ), (48a) 

Gi cote ldGi ,, fi ,v 

Si = 48b 

r r at) 

where i = P, H, T, and R. Eq. (H7|) has a specific solution: 



x- 



r 



Ui = \ I x 2 Sidr. (49) 



o 



3.3. Uniform Rotation Rate 

3.3.1. Uniform rotation equipotential surface 

We want to use this special case to show how to obtain the equipotential surface x = 
x(m, 6). 

For rotation at the angular velocity Qz, we can use Eq. (|49p to calculate the radial 
gravitational acceleration perturbation Ur. The result is 

Gr = Tle = ~0 2 r sin 29, 
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Sr 
U r 



3^2 



1 



2 fi 2 (cos2fl + -), 

--^ 2 r(cos20 + -) 
2 v 3 J 



Here we assume that Q = 0(r) does not depend upon 9. Eq. fj4*0l shows that we need 

U R -Tl r = -tt 2 r(cos29 + l). 

As the first approximation, we assume p/p m — 1, x — 1, and Up = Ud = in Eq. (j4"Ul) . For 
a slow rotation in the sense that the centrifugal acceleration Q 2 r e is much smaller than the 
corresponding gravitational acceleration Gm/r 2 , we obtain 



x 



(0) 



1 - -a (cos26> + 1), 



(50) 



where 



n 2 r 3 

Gm 



a = 



We can further improve the result by taking into account w = p/p m , A = (%-) , ?7p, 
and % in Eq. (SO]): 

1, 



A (o 



( O 



5 



(o 



5 



1 - ia (cos26> + 

^a sin 26* , 

1 Gma . 
-- — —sin 20, 
2 

1/ oflJ l,3Gma 
2 (C0S ^ + 3 ) ^^' 

-i(cos20 + ^)47rGa o p, 



cos20 + i)^ O) , 



where 



l(0) 



-~(cos20 + ~)&g\ 



^rfr, bp* = AnG I a pdr. 



The corrected equipotential surface function is 



x 



(i) 



1 + l[ao + ^(&g } - &£>)] - ^(008 20 + r 



(51) 
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where 



1 T 

ai = a + a' , a' = -[a + Q^( b D b p)\- 



According to the definition of x, it should equal unity at the equator. This requirement fixes 
the expression of x as follows: 



x 



(i) 



1 - ^-a 1 (cos2# + 1). 



(52) 



From now on, we shall show this form only, which will be referred to as the normalized form. 
Using Eq. (1521) or its non- normalized form we can improve w, A, Up, and Up,: 



s 



s 



U 



U 



l-^a 1 (cos20 + i; 



— ai sin 29, 
2 

iGmax . nQ 
- — sin 20, 

2 r 2 

1 . „ 1 . 3Gmai 
-(cos2* + -)— ^, 

-^(cos20 + -)^Ga lP , 

~(cos20 + ~)&« 

-h CO s26 + l)b%, 



where 



■ [ — - dr, bp* = AtcG f a x pdr. 
Jo r Jo 



The more accurate equipotential surface is thus expressed by 



x 



(2) 



1 - -a 2 (cos2# + 1) 



where 



1 r 2 

a 2 = a + a[, a[ = -[a x + ^— ( 6 2 ~ b p')\- 



(53) 



To keep iterating, we find the following recurrence relation for i = 1, 2, 3, 



x 



1 - -ai(cos26» + 1), 
4 



(54) 
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(55a) 
(55b) 

(55c) 
(55d) 

Using the equipotential surface profile, Eq. fl5^|) . we can calculate the following quantities: 



Q {i) 


= 4vrr e 2 p(l - ^Oi), 


(56a) 




= 1 - ~a;(cos20 + -), 


(56b) 




= -ajsin26 l , 


(56c) 




1 Gmai . 
= 2 r 2 sm2£? > 


(56d) 




= i(cos2^ + i)6«, 


(56e) 


rr(«) 


= -^(co S 2e + i)&«. 


(56f) 



The gravitational acceleration perturbations due to rotation are 

= ~(OV-^)sin20, (57a) 
2 r l 

U {i) = -J[3Q 2 r + (&g } -6g)](cos 2 + J). (57b) 

With inclusion of the rotation effects, the gravitational acceleration vector can be expressed 
as follows: 

gf = ^-\[3Q 2 r + (b i S-bf)](cos2e+ 1 -), (58a) 
gf = i (n 2 r~^p^ sin 26. (58b) 

We know r = r e x^ . Since b^p and bp are integrals over r from to r, we know that the 
gravitational acceleration perturbations £/w and do not vanish outside the star. 



where 



di = a + a-_ l5 
6 « = 3 f r Grna* 



'o 

rr 

b { £ = AttG / a iP dr. 



3.3.2. Uniform rotation-like magnetic equipotential surface 



Rotation has a global velocity field v = (0, 0, Qr sin#). We can choose a toroidal mag- 
netic field B = (0, 0, (inpY^Qr sin 9) to mimic rotation at the rate f2. We use this magnetic 
configuration to show the calculation method for the magnetic equipotential surface and the 
difference between rotation and magnetic effects. 

The first step is to calculate two components of H: H r and Tig. They are 

H T = -tt 2 rsm 2 9, 
Hg = —Q 2 rsm9cos9. 

Comparing them with the corresponding TZ r and 1Z$, we can see that their signs are opposite. 

We also need the plasma (3 parameter. Its definition is the ratio of the total pressure 
Pt over the magnetic pressure P m = |pf2 2 r 2 sin 2 9. Using j3 = 2P T / p Q 2 r 2 , we have (3 = 
(3q/ sin 2 9. Magnetic pressure causes a density change. The density (p/po) with/without the 
magnetic field is related to each other by the formula p = po/(l + 1//3), or p = po/(l + 
c 2 sin 2 #), where we have used c 2 = l//3 to replace /3 . We know c 2 = por 2 Q 2 /2P T . 

The next step is to use Q H = Hg to obtain the source term Sh' 

S H = ^n 2 (cos29 + ^). 

Substituting it into Eq. f|49|) . we obtain 

U H = ^Q 2 r(cos29 + ^). 

Using Ujj and 7i r in Eq. (|40p . we obtain the first approximation to the magnetic equipotential 
surface 

= l + -a (cos2fl + l). (59) 

Comparing Eqs. (1301) and (!59j) . we can see that the oblateness e = (r e — r p )/r e = ±a/2 is 
positive for rotation, but negative for magnetic fields, where r p is the polar radius. 

The following steps differ from the rotation case since the magnetic effect on density, 
which comes from the integral p m , cuts in. The density correction to the equipotential surface 
can be expressed by c 2 in the recurrence relation 



x« = l + -a i (cos2fl + l), 



(60) 
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where 



«0 + -C2 + CJj 



i-l> 



6 (i> 

p 



l(0 



Gma, 



47rG / (a» + c 2 )podr. 

'0 



(61) 
(62) 

(63) 
(64) 



Using the equipotential surface profile, Eq. fl60|) . we can calculate the following quanti- 



ties: 



(<) 



17 



1 2 

4vrr 2 p (l + - -c 2 ), 

1 + -(a i + C 2 )(cOS20 + ^ 

— aj sin 26, 
2 

— sin2f, 

2 r 2 

1, _„ 1 



(i) 



(cos20 + -)&£ 



U 



(cos 20 + ^)&g ) . 



1. 

3' 



(65a) 
(65b) 
(65c) 
(65d) 
(65e) 
(65f) 



Since the magnetic effect on density has been totally absorbed into c 2 , the integrant in the 
integral b£ involves po, instead of p = po/(l + c^sin 2 ^), which is the same as above. The 
gravitational acceleration perturbations due to a rotation-like magnetic field are 



g(i) 



1 ( ^ Gmai\ . 
-- (fi 2 r--^j sin 20, 

~(3QV + &g ) -&g ) )(cos20 + ~). 



(66a) 
(66b) 



Including the rotation-like magnetic effects, we obtain the expression for the gravitational 
acceleration vector: 



9 



Gm 1 
2 



2 + ^(3fi 2 r + &«-&g ) )(cos20 + i), 



9e 



(i) 



1/^9 Gma; . 
— ft 2 r — 1 sin 20. 



(67a) 
(67b) 
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3.3.3. Uniform rotation-like turbulent equipotential surface 

Solar turbulent data are given by the three-dimensional (3D) numerical simulations 
within a small volume that contains the super-adiabatic layer (SAL) of the Sun. The turbu- 
lent pressure P t = -^pv"v" peaks at the peak of SAL. The peak value is about 17% (Robinson 
et al 2003; Stein & Nordlund 1998). Since the simulations are restricted to a small range of 
the colatitudinal coordinate and all the turbulent velocity components are the averaged ve- 
locity variance over the colatitudinal coordinate, the ^-dependence of the turbulent velocity 
is unknown. Turbulent velocity may have two components, one is ^-independent, and the 
other is ^-dependent. The latter must be much smaller than the former. 

The ^-independent component has nothing to do with the equipotential surface, but the 
^-dependent component affects the equipotential surface. In order to address the difference 
among rotation, magnetic, and turbulent effects, we assume that the ^-dependent component 
of v"v" equals |f2 2 r 2 sin 2 9, and that of v'gVg equals zero or Q 2 r 2 sin 2 9. As a result, we have 

% = Ttt 2 rsm 2 9, (68a) 
% = ±Q 2 rsin9cos9. (68b) 

It is interesting to note that the signs of both lZ r and TZg are the same ("+"), those of 
both Ti r and Tig are the same ("-"), but those of % and T e are opposite to each other ("=f" 
vs "±"). We have shown above that the sign determines the sign of the oblateness of the 
equipotential surface. We thus anticipate something new for turbulence. Following the same 
procedure as obtaining Eq. (160]) . we obtain 

x (0) = 1 ± i(2a )(cos2fl + 1). (69) 

The new outcome is that the coefficient doubles, here = Q 2 r 3 /Gm as above. The recur- 
rence relation thus becomes 

x (i) = l±iai(cos2# + l), a i = 2a ±^c 2 + a' i _ 1 , (70) 

where ft = l/c 2 sin 2 # is the turbulent ft parameter. The expression for a[ is the same as 
above. 

When we assume that the ^-dependent component of v'gVg equals twice that of v"v", 
we obtain the same gravitational acceleration as that for rotation, Eqs. ( 158aP -( l58b|) . except 
that is defined in §3.3.21 when we assume that the ^-dependent component of v'gVg equals 
zero, we obtain the same result as that for the rotation-like magnetic field, Eqs. (I67ap - fl67bj) . 
Therefore, turbulence plays a role of either rotation or magnetism. The criterion is: we have 
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the rotation/magnetism effect when the transverse turbulent velocity is larger /smaller than 
the radial turbulent velocity. 

Solar 3D turbulence simulations show that the transverse turbulent velocity is smaller 
than the radial turbulent velocity near the solar surface. We thus expect some magnetic 
effects therein. 



3.3.4- Uniform rotation-magnetism-turbulence equipotential surface 

In the general case, we can express the equipotential surface in the same formula as the 
magnetic equipotential surface: 



where 



an 
aT 
a H 

a- 

Up 
°D 



x {i) = 1 + -a;(cos2fl + 1) 



1 



a H ± 2a T - a R + ~(c m + c T 2) + a-_ 1; 

Gm ' 
Gm ' 
Gm ' 

Gmai 
— „ — C", 



47rG / (a* + c^ 2 + c T 2)Podr. 
Jo 



(71) 

(72a) 
(72b) 
(72c) 
(72d) 

(72e) 
(72f) 

(72g) 



Using the equipotential surface profile, Eq. (loTij) . we can calculate the following quanti- 



ties: 



(i) 



47rr e p 



1 2 . . 

1 + -Oj - g( c ^2 + CT2J 



— 



1 + + c ^2 + c T2 )(cos26' + -), 
— —at sin 2$, 



(73a) 

(73b) 
(73c) 
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uf = 



lGmai . 
— sin 26', 

1 



2 r 
1 



U { S = ±(cos26> + ±)&g>. 



(73d) 
(73e) 
(73f) 



The gravitational acceleration perturbations due to rotation, rotation-like magnetic field 
and turbulence are 



0(0 



(n 2 H -n 2 R ± n 2 T y - 



Gmaj 



sin 20, 



l[3(Q 2 H -n 2 R ± n 2 T )r + b { $ - &?](cos20 + \). 



The gravitational acceleration vector in the system is: 



W = ^ + lm 2 H -n 2 R ±n 2 T )r + b^-b^}(cos29 + l), 



9r 



r 



1 Gmai 



2 r 2 



sin 20. 



(74a) 
(74b) 

(75a) 
(75b) 



So far we have assumed that fli (i = R, H, T) are uniform. They depend upon r and 9 in 
general. This is so-called differential rotation. We deal with the more complicated situation 
in the next section. 



3.4. Differential Rotation Rate 

3.4-1- Differential rotation equipotential surface 

Not all form of differential rotation is non-singular. Whether some differential rotation 
is singular is determined by Sp, which contains the term Qp cot 9. This term is non-singular if 
Qp has a sine function factor, sin 9. This criterion yields the following non-singular differential 
rotation profile: 

N 

Q 2 (r, 9) = J2 Q 2n(r) cos 2n9, (76) 

n=0 

where N is an finite integer. This form of expression for Q 2 is physical because physical 
solutions should not be singular. 
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The first order of approximation to the equipotential surface is 

1 -iE a 2°n ) [ cos2 ^ + (- 1 : 



a; 



(0) 



\n-ll 



(77) 



n=l 



where 



,(o) 



1 r 3 



(0) 
2)i 



4 Gm 
1 r 3 



[2(3fi - fio) + 2(^2 + fi 2 ) " («4 + ^)], 

{[(2n + l)n 2n _ 2 - n 2n _ 2 ] + 2(H 2n + fi 2n ) - [(2n - l)n 2n+2 + fi 2n+2 ]}, 



4 Gm 

u R = -^{(2n +n 2 ) + (6n + 2n 2 -n 4 )cos2^ 

+ [O + l)fia«-2 + 2H 2n - (2n - l)n 2n+2 ] cos2n6}. 

n=2 

We have defined f2 = 7 /J" ^W r ; etc. 

The next step is to calculate Up^ and which are used in Eq. ([30]) . They are 



CE7 



(0) 



l (0) 
'-2n 



n=l 



cos 2n# + 



(2n - l)(2ra + 1) 



AT+l 



A (0) 



\ 5^[n4° ) sin2n0], 



71=1 



AT+l 



C7 



(0) 



1 Gm \ . (0) . „, 

71=1 

^ AT+l 

2 5Z 6 pL cos2nfl, 

n=0 
AT+l 



n=0 



where 



l(0) 
7 P0 



-(0) 
^P2n 



Gm 



N+l 



11 ' „-!. 

Gm 



^na^dr, 



N+l 



n(2n + l)4 n } + Y 2ka 



(0) 

2A: 



fc=n+l 
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6 (o) 



N+l 



i {0) 



o r f^i (2n-l)(2n + l) 



dr. 



41 = 4ttG / pagdr. 



The corrected equipotential surface function is 

JV+l 

4 



x 



(0 



1 AT+1 

i- I E a £[cos2^ + (-ir i ] 



n=l 



where 



,(0 



.0' 



Q w = 4vrr^ 



(0) (i-l)' 



^ 5 



1 \ - c 

1 + - > 7 

2 ^-f (2n - 1 



.(0 

2n 



1 iV+l 



(7* 



(79a) 
(79b) 

(79c) 



-(2n-l)(2n + l) 2^ 

for £ = 2, 4, 6, ■ • 2(N+1), and i = 1, 2, 3, - - • 

Those terms with I ^ 2 in Eq. flTHj) are pure differential rotation effects. The term with 
t = 2 also contains some differential rotation correction. 



3.4-2. Differential rotation-like magnetic equipotential surface 
The following toroidal magnetic field mimics the differential rotation, Eq. fl76l) : 

Bt(r,9) = (An p) 1/2 n(r,6)r sin 6. 
The system has the following equipotential surface: 



(80) 



x 







1 N+l 

l + iE a £[cos2^+(-ir i ] 



5i) 



where 



,0 



,0' 



n=l 



) + -q 4 " 
2l a £ Gm °pi)\i 



(i-iy 

i i 



(82a) 
(82b) 
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for t — 2, 4, 6, • • •, 2(N+1), and i = 1, 2, 3, • • •. The starting point af^ is the same as 
above except that Uh = —Ur. The coefficients q are defined by the relation p = po/(l — 
\ YZ=l °2n sin 2n9). They are 



c 
c 2 
ct 



p r 2 1 



2P T 2 
p r 2 1 



-(2^ -0 2 ), 



2P T 2 
p r 2 1 



2P T 2 



-(2^0-2^2 + ^4), 

(£lg-2 — 2Vle + 



for i = 4, 6, 8, - 2(N+1). 

Using these expressions, we can calculate the following quantities: 



Q w = 4nr 2 e p 



JV+1 



1 1^ 



(») , 
°2n + C 2n 



JV+1 



2 2 ^ (2ra- l)(2ra+ 1) 2 ^ 

n=l K ' v ' n=l 



«2r, 



JV+1 



(0 



n=l 
JV+1 



cos 2n9 + 



(2n-l)(2ra+l) 



A« = 



- ^[nag sin 2n0], 



2 

1 Gm 



n=l 
JV+1 



J^[na£ sin 2n0], 



17 



(0 



17 



(0 
jj 



2 r 2 

n=l 
JV+1 

Jp2n cos 2n#, 

n=0 
JV+1 

2 



- E>?> 



-, JV+1 

9 u D2n 



COS 2n#. 



n=0 



(83a) 
(83b) 
(83c) 



The coefficients are the same as above, but coefficients bp are different from above. They 



are: 



-'do 



&« = 47rG / po ( a w + Q ) rfr for £ = 2, 4, 6, • • •, 2(N+1) 



„ r /V+l , (j) X , 

7o 2^(2n + l)(2n-l)' 



,(') 
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3.4-3. Differential rotation-like turbulent equipotential surface 



The differential rotation-like turbulent parameter is the same as Eq. f[76|) . This system 
has the following equipotential surface in the first approximation: 



JV+l 



X 



(0) 



lT~X>£![cos2n0 + (-l) 



n+l] 



n=l 



where 



4°> 



L 2n 



1 

4 

1 r 3 



[2{3Th + o ) + 2(n^ - n 2 ) - (o 4 - n 4 )], 



4Gm 
1 r 3 



4 Gm 



{[(2n + l)^2n-2 + ^2n-2) + 2(fi 2n - ft 2 n) - [(2/1 - l)0 2n+2 - fi 2n+2 ]}, 



C/ T = T -{(2n + fi 2 ) + (6fi + 2fi 2 -fi 4 )cos2^ 



N+l 



+ [(2n + l)ft 2 „_ 2 + 2fi 2n - (2n - l)fi 2n+2 ] cos2nfl}. 



n=2 



We have the following recurrence relation: 

. JV+l 

1tT«SJco S 2^ + (1. 



(0 



\n+ll 



55) 



where 



(0 



n=l 



,(0) 



a 



(86a) 
(86b) 



+ 

for £ = 2, 4, 6, • ■ -, 2 (N+l), and i = 1, 2, 3, ■ • -. We can use it to express the following 
quantities: 

JV+l 



Q<" = 4,wr 2 ePo 



1 + 5»±5E 



AT+1 (i) 

a 2n + C 2n 



(0 



2 ^-f (2n- l)(2n+ 1) 2^ 

n=l n=l 

l + -^(a£ + c 2 „) cos2n6+ ' 



2V+1 



n=l 



(2n - l)(2n+ 1) 



n=l 
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U® = 



u 



(0 

D 



1 n N+i 

n=l 

n=0 
N+l 



4e» 



2 J -D2ra 
n=0 



cos 2n#. 



5.^.^. Differential rotation-magnetism-turbulence equipotential surface 



Put all three sources together, we have the following recurrence relation: 



N+l 



x® = l + \ y Ea®[<xx2nB + (-in 



(87) 



71=1 



where 



a 



(0 



in „<°> -U " /V __l_ ^ ^ _l_ K' 

a H£ ' a Tl ~ a M + 2 £ 

„2 



(i-iy 



4 = 2^ +^( 6 w-y> 



(88a) 
(88b) 



for £ = 2, 4, 6, • • •, 2(N+1), and i = 1, 2, 3, • • •. The coefficients b® are the same as above, 



but coefficients b® are: 



JV+l / (i) 



Po(4n + C ^2n + C T2n )dr 



, — (2n+l)(2n-l) ' 

4ttG / p (a? + c m + c Te )dr for £ = 2, 4, 6, • • •, 2(N+1). 
Jo 



The useful quantities are: 



Q(0 = 47rr e 2 po 



1 1/ , \ 1 a 2n + C #2n + C T 2n 1 / 

1 + 2^o + c T0 ) - 2 L (2n-l)(2n+l) " 2 ° 2 " 

n=l v 7 v ' n=l 



CE7 



j N+l 



n=l 
JV+1 



cos 2n6 + 



(2n- l)(2n + 1) 



A (i) = -J5>42sin2n0], 



n=l 
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-I n N+l 

n=l 
1 N+l 

n=0 



2 

n=0 

Cm 



JV+l 



n=0 

N+l 



gf = -^(^ T ^-^)rsin2^ + ^^na«sin2^. 

n=l 

This is the general case for rotation, the rotation-like toroidal magnetic field and turbu- 
lence. The recurrence relations given here reflect the real cause-effect relation. The source 
terms {Ur — TZ t ), {Uh — 'Ht) an d {Ut — %) are the causes, and Up, Uo and w are their effects. 
When some asphericity sources are present, the spherically-symmetric star should readjust 
to assume an aspherical equilibrium configuration. The recurrence relations describe the 
readjustment procedure. 



4. METHOD OF SOLUTION 

4.1. 2D Stellar Structure Equations with an Known Equipotential Surface 

For the cases studied above, we can use the recurrence relations to calculate the equipo- 
tential surface functions to certain accuracy. The result is denoted 

(°°)_ From 

now on, we use the un-superscripted symbols to express the corresponding limits, for ex- 
ample, ag = a^, and so on. We then use x to calculate functions w, A, Q, etc. This is 
equivalent to solving the Poisson equation for the gravitational acceleration vector. 

With the help of the equipotential surface, what we need to numerically solve for are 
r e , Pt, T, and L, which are governed by the following four equations: 

£ = jf, (89a) 
os Qr e 

dP' Gm 2 . , 

aT = (89b > 
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ar 


dp' 


ds 


ds 


dL 


mm 


ds 





V ra d radiative 
V r convective 



dt J L Q r e p 
Here r' e = lnr e , r = r e x, w = p/p m , A = (dx/d9) m , and 



1 pCpTl m v conv Vg 



ad 



2 1 + v conv / v r 



39c) 



T dSr\ _ (ggd) 



dF 

* = F e cot# + -/, (90a) 

00 

F e = P(F l + F 2 + F 3 ), (90b) 
P = GmQxA (90c) 

= (90d) 

^ 2 _ 1 pC P Tl m v conv V 
2 1 + v conv /v r ' 



(90f) 



The variable \I/ has a term that is proportional to the following expression: 



x <9 A 

A = A cot 6 + — 

N+l 



N+l 



- na ?n 



71=1 



n=l 



JV+1 



n(2n + l)a 2 „ + 2/ca 2 fc 



k=n+l 



cos 2n#. 



(91) 



The second term is FF 2 , where F is defined in §A.ll The required \I> is the sum of these 
two terms: 

GmQA x 2 3 ~ 2 



(F + F + F ) + FF 



47rr3F T 

The other supplement quantities are given in §3.4.41 



(92) 



4.2. Linearization of 2D Stellar Structure Equations 

The construction of a two-dimensional stellar model begins by dividing the star into M 
mass shells and N angular zones. The mass shells are assigned a value Sj = logmj, where m; 
is the interior mass at the midpoint of shell i. The angular zones are assigned a value 8j. A 
starting (or previous in evolutionary time) model is supplied with a run of (F/, T[-, r[, Lij, 
U i:j = 0, Qij = 0) for i=l to M and j=l to N. 
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Different terms in Eqs. (I89ap - fl89dj) have different derivatives with respect to the stellar 
parameters (Pt, T, r, L). These derivatives are needed to write down the linearized difference 
equations. We hence rewrite them as follows: 



dP' 

ds 
dT 

ds 
dr' 

ds 

dL 

ds 

The symbols used above are defined as follows: 



V, (93a) 

T, (93b) 

K, (93c) 

IX ( 93d ) 



Gm 2 , . 

P = 94a 

T = VV, (94b) 

ft = (94c) 

, raw ( ^c/SVA 

£ = -r- e - T -r- > ( 94d ) 



L & \ dt 
raw 



C 2 = --^-F e cote, (94e) 
L & rp 

< 3 - (94f) 



We use the central difference scheme to approximate the stellar structure equations. 
The corresponding difference equations are 

Pp = (P;-PU)-^As l (V l + V^ 1 )=0, (95a) 
F% = (Tr-TU J )- 1 -As t (T lJ+ T l _ lj ) = 0, (95b) 

Fr = {m i -R' i _ 1 )-^As i {n i + n^ l ) = ^ (95c) 

i 3 

F$ = (L y -L i _ li )--As i 53(4 + £{_ li ) = 0, (95d) 

for i = 2 to M, and j = 1 to N. The linearization of Eqs. (I95al) - fl95dl) with respect to (SP(j, 
ST^, 5r' ip and SLq) yields 2(M - 1)N + 2(M - 1) equations for the 2MN + 2M unknowns. 
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The N + 1 additional equations are supplied by the boundary conditions at the center: 

Fr = r[ - [ Sl - ln(Q/3)] = 0, (96a) 

Fl j = ^,-^4=0, (96b) 
e=i 

where j = 1 to N. Another N + 1 additional equations are supplied by the boundary 
conditions at the surface: 

Fr +1 = R'm ~ a \P'u ~ a 2T' MN — a 3 = 0, (97a) 
pM+ij _ L> M .(i n L M ._ a4 p^_ a5 T^.-a 6 ) = 0, (97b) 



where j = 1 to N. The F equations are linearized, 

/ fj jpij f) rpij fj rpij fj rpij 

- EE 4«i^i.^^4«i 

1=1 k=l x 1 1 lK 

for w = T, L, (98a) 
= E ( s]f 5 ^ + 9]f 6P I) for w = R, P , (98b) 



i=i 

where % — 1 to M; and j = 1 to N. The summation over 1 has non-zero terms only for 1 
= i-1, i; the summation over k has non-zero terms only for k = j. See appendix A for the 
coefficient matrix elements. 

Since we explicitly take advantage of the equipotential surface function x, we can express 
the derivatives of all dependent variables with respect to 9 in terms of A, which is the 9- 
derivative of x on the equipotential surface. This unchains the explicit binding between 
adjacent angular zones and allows us to treat each zone as if it is a one-dimensional problem. 
However, the implicit binding cannot be broken because of the mass conservation requirement 
that is characterized by the parameter Q, which is an integral over all zones. 

These equations can be solved by means of the Henyey method. 



4.3. Non-equator Reference Surface 

So far we have used the equator as the reference surface. This is not necessary. We can 
use the other reference surface instead, say, 9 = 9 . The equator is only a specific example 
where 9o = ir/2. We need a non-equator reference surface when the applied field peaks at or 
near the equator. We use the subscript "f" as the indicator of the reference surface 9 = 9 . 
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Since r = rfX, the equipotential surface x = x(rf,9) should be normalized to unity at the 
reference surface 9 = 8 . We give different formulas as follows: 

, N+l 

x {i) = 1 + - [cos 2n9- cos 2n9 ], (99a) 

71=1 

f 1 1 AF+1 

L n=l 

We use subscript "f to replace "e" in the other formulas and/or equations. 

5. HIGH-PRECISION 2D SOLAR MODELS 

The solar variability models need to be accurate enough to match the seismic structures 
of the Sun (Gough et al 1996), as the (ID) standard solar models do (Bahcall et al 2006 
and references cited therein). Standard solar models (1) use the most accurate available 
input parameters, including radiative opacity, equation of state, and nuclear cross sections, 
(2) include element diffusion, and (3) have a high numerical resolution. Our 2D models 
inherit all these features because our 2D code described in this paper is a natural extension 
of YREC (Yale Rotation Evolution Code) to two dimensions. We also tested its ID counter- 
part with turbulence (Li et al 2002) and made sure that the resultant ID solar models are 
accurate enough to meet with our accuracy requirements. We further tested the ID code 
with magnetic fields and turbulence (Li et al 2003) to make sure that it is accurate enough 
to discern the solar cycle-related p-mode frequency changes. These demonstrate that the 
first dimension is accurate enough to discern the solar cycle-related changes. The number of 
mass layers used in both ID and 2D model calculations is more than 2500. 

5.1. Error Controls 

Here we describe how we control the numerical errors to meet with our accuracy re- 
quirement. 

5.1.1. Radial 

This is the same as its ID counterpart. The numerical errors are controlled in terms of 
two parameters e F and ec' 

\Fij\ < e F , and \5w ij \ < e c (100) 



(i) 
l 2n 



CH2r 



CT2n 



(2ra- l)(2n+l) 



+ a^l cos 2n9 



.(99b) 
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for i = 1 to M+l, j = 1 to N, and w = P, T, R, and L. See Eqs. 095ati97bl) for the definition 
of J*?- 

The ID standard solar models have ep ~ ~ 10 -6 , which is the relative accuracy of 
the numerical solution of the stellar structure equations. We use the same values of e for our 
2D solar models. 



5.1.2. Colatitudinal 

From §4. II we can see that the colatitudinal factors affect the stellar structure equations 
in terms of x, Q, w, A, and A. The quantities Q and w are the integrals of x over 9, and 
A and A are the (first-order and second-order) derivatives of x with respect to 9. Therefore, 
the colatitudinal errors are determined by the error of the equipotential surface function x, 
which is defined by Eq. (1401) . 

For rotation, rotation-like magnetic field, and/or rotation-like turbulence that are sym- 
metric with respect to the equator, Eq. (j4"0~|) can be rewritten as follows: 



/ 1 N+l \ -V 2 / oo \ V 2 

x = q- 1 ' 2 I 1 - -^c 2n cos2n0 J 1 + x 3 ^ u 2n cos 2n8\ . (101) 

\ n=0 / \ n=0 / 

In doing so we have rewritten p m (Eq. HTJ and p as follows: 



Pm = Pox 2 q, 



( 1 N+1 

p = p ( 1 - - ^ C 2n COS 2n ® 



n=0 



We have also used the fact that U oc x, Ti r oc x, % oc x, and 7£ r oc x. The quantities used 
here are defined as follows: 



g = / j 1 — c 2 „ cos 2n9 J x 2 sin 0cf0, 

«2n = 4n + oTTM^Sn - 6p2n), (102b) 

2 Gm 

6.0 = iTrG^ ^ (2n + 1)(2n _ ir (102c) 

= 4vrG / p (fl2n + c 2n )dr, (102d) 
Jo 



'102a) 



6 
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?P0 



J P2n 



Gm \ -» 

— — na 2n dr, ( 102e) 

1 T * n=l 



Gm 



o ' e 



n(2n + l)a 2n + 2/ca 2fc 

fc=n+l 



dr, (102f) 



for n = 1 to oo. Here we have used the Fourier series to express the normalized equipotential 
surface function x: 

oo 

a; = l + ^a 2n [cos2^ + (-l) n+1 ]. (103) 

71=1 

For pure rotation, c 2n = for all n (n = to oo). 

In practice, we have to truncate the infinite Fourier series to approximate x, 

Af 

XN = l + ^[cos2n#+ (-l) n+1 ]. (104) 

71=1 

Since | cos2n#| < 1, the truncation error can be estimated as follows: 

oo 

e x = \x - xm\ < 2j l fl 2n|- (105) 

n=N+l 

If the a 2n 's are rapidly decreasing, which is the typical case, then the truncation error is 
dominated by g^jv+i)- We can thus use a 2 (N+i) as an estimate of the truncation error of x: 

£x ~ 1 02(^+1)1- (106) 

We want to achieve a relative accuracy of 10~ 6 for the stellar parameters P, T, R and 
L in the 2D model, the same as in the ID standard solar model. This requires the similar 
relative accuracy for x. Since x is of the order of magnitude of unity, its relative error is 
the same as its absolute error. In order to achieve such high an accuracy, we use three-level 
iterations to solve Eqs. (I101H103P 

The first-level iteration is given in §0 in terms of the recurrence relations, which are 
based on the linear approximation of Eq. fj40l . The convergence criterion is |a 2 „ — a 2n ^| < e 
for i = 1 to N+l, where e = 10 -6 . The converged denoted by a 2n . The second- and 

third-level iterations are used to do nonlinear corrections. 

The second-level iteration uses 

Af+l 

x f) = l + J2 4n[cos2n9 + (-l) n+1 ] (107) 

71=1 
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as the initial guess for x in Eq. (110 II) . The updated xfj is normalized as follows: 

x fj = x fj - xfj (6 = tt/2) + 1 (108) 
for i = 1, 2, 3, The convergence criterion is 

\4i-4r 1] \ <e- (109) 

The converged xfj is denoted by xn, which is then expanded as the Fourier series to prepare 
for the third-level iteration: 

oo 

x n = 5^a^ COS 2710. (110) 

n=0 

We have to truncate Eq. (1 11 01) to go further. The truncation criterion is 

|«2Ar| > e and \a 2n \ < e for n > JV + 1. (Ill) 

Generally speaking, jV > JV + 1. 

Using (n = 1 to A/") as the initial guess for a^ 1 , denoted as bfj, we repeat the 
second-level iteration to update b 2 l- The convergence criterion is 

\bfl-bt l) \<e (112) 

for n = 1 to M . The converged bfl's are denoted as a^ 1 ■ Using a^ 1 , we can calculate x, Q, 
zu, A, A, and other quantities such as g r and gg. 

Extensive numerical experiments reveal that the dominant error sources come from 
Eq. (1102f[) . whose integrand is proportional to the Fourier expansion coefficients of A, 
Eq. ([91]): 

Af+l 

A 2n = n(2n + l)a 2n + V] 2ka 2k for n = 1 to Af+1. (113) 

k=n+l 

Its first term originates from the second derivative of the equipotential surface x. The 
coefficient n(2n + 1) of a2„ in the first term will substantially magnify the error of a 2n when 
n is big. In order to control this error, we calculate the maximal value of the ratio of the 
centrifugal over the gravitational acceleration for pure rotation, denoted as rj, we define r\ 
as the maximal value of 1/(3 for magnetic fields and/or turbulence. Numerical experiments 
show that the convergence criterion is e = max(ei?, ec, if). 
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5.2. Examples 

5.2.1. Uniform rotation 

This is the simplest case. First of all we calculate a high-precision (ID) standard solar 
model by using the convergence criterion = ec = 1 x 10~ 10 . We use it as the benchmark. 
We then use zero-rotation rate (Q = 0) to calculate a series of 2D solar models by using the 
convergence criterion e = €f = €c from 1 x 10~ 3 to 1 x 10~ 9 . The numerical accuracy of the 
2D solar models is measured in terms of their relative errors with respect to the standard solar 
model. The model is represented in terms of runs of pressure, P = P(m,9), temperature 
T = T(m,9), radius r = r(m,9), luminosity L = L(m,9), and density p = p(m,9). The 
numerical accuracy of the 2D solar models is thus defined as the maximal value of the relative 
errors for all five variables over all grid points. The results are shown in Fig. [TJ in which the 
symbols mark the data points. The figure shows that we can achieve a precision significantly 
better than 1 x 10~ 6 , which is accurate enough for the relevant solar applications. Since we 
avoid numerical derivatives and integrals, the results are independent of the grid size in the 
second coordinate 9. This is confirmed by the detailed model calculations by setting N = 9, 
17, and 33, where N is the number of grid points in the second dimension. For both ID and 
2D models the first dimension has the same grid point number M = 2576. 

When the rotation rate is nonzero, i.e., Q ^ 0, the relative differences between the 2D 
and ID models such as Sp = [P(m, 9)—P(m)]/P(m) etc can be considered to be the rotation 
effects. They are functions of the rotation rate Q, convergence criterion e, the mass coordinate 
m and colatitude coordinate 9, for example, Sp = Sp(m,9;Q,e), Sp = ST(rn,9;Q,e), and 
similar expressions for r, L and p. Their accuracy is estimated by the corresponding value at 
the zero-rotation rate. Fig. [2] shows how the maximal value of Sp, Sp, S r , Sp, and £ p changes 
with Q, where we fix e = 1 x 10~ 6 (solid line) or 1 x (dotted line). So the relative error 
is of the same order indicated by the dashed line (e = 1 x 10 6 ) and the dot-dashed 

line (e = 1 x 10~ 7 ) in the figure. 

To see where the maximal rotation effect takes place, we plot £r = £(m(R), {9}; Q) as 
a function of R/R Q and Q in Fig. [3], where £p is the maximal value among £ P , £ T , £ r , £ L 
and £ p over all zones, and R is the radius of the mass shell m in the standard solar model. 
Similarly, we have S$ = £({m}, 9; Q). Since it changes little with 9, we do not need to plot 
it. Fig. [3] shows that the maximum takes place at the base of the convection zone or near 
the surface. Fig. H] shows the detail dependence of £p, £t, £ r , £l and £ p on R/R Q and 9. It 
also shows the equipotential surface x, Fg, 5g r and gg, which have no ID counterparts. 
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5.2.2. Uniform rotation-like magnetic field 

The uniform rotation-like toroidal mag netic field is B = (0, 0, (Anp) 1 ' 2 ^ sin 6>) . We 
repeat the similar model calculations to rotation. Figs. [5] -13 show the results. Once again, the 
high-precision is achieved. Comparing them with Figsj2}H] we can see rotation-like magnetic 
fields affect stellar structures in a different way from the rotation: magnetic effects take place 
in the convection zone and peak near the surface. Rotation-like turbulence behaves like a 
rotation-like magnetic field. 

5.2.3. Differential rotation-like magnetic field: torus 

The torus field is a rotation-like toroidal magnetic field, B = (0, 0, (AnpY^VLr sin 9). 
The magnetic rotation rate Q is defined in Appendix [Bj There are two torus tubes that are 
parallel to the equatorial plane since they are assumed to be symmetric with respect to the 
equatorial plane. As a result, there are four circles on any meridional plane. 

Unlike the uniform rotation rate, we should first find out the discrete Fourier transform 
of the square of the differential rotation rate Q, Q 2 , which is equally discretized in the range 
of 9 from to ir/2, namely Q{ for i = to N. Here N should be a power of 2. We calculate 
O 2 in the first quadrant and then extend it to the other three quadrants according to the 
symmetry described above. Its discrete Fourier transform F n are finally calculated by means 
of the Fast Fourier Transform (FFT) of a real function (See the subroutine realft.for given 
in Numerical Recipe) for n = to 4N. Each pair of the data contain the real and imaginary 
parts of the FFT except for the first pair. The imaginary part vanishes since Q 2 is a real 
function of 9, which is now in the range of to 2tt. The odd components vanish due the 
equatorial symmetry. We use y n to denote the nonzero components. The nonzero F n contains 
Fq, which is twice the uniform component, yo = F /2; and F\, which stores the twice of the 
Nyquist critical wavenumber component, y^ = Fi/2; and the even components y n = F 4n for 
n = f to N-f . Consequently, we have 

N 

Q 2 = ^2y n cos2n9. (114) 

n=0 

Fig. [S] contains nine sub- figures for the Gaussian profile defined in Appendix B, in which 
fio = 3x 1CT 5 . Sub-figure (1,1) shows the reciprocal of the plasma (3 parameter as a function 
of (R/Rq,9), which is defined as the ratio of the gas pressure over the magnetic pressure: 
1/(3 = \pVL 2 r 2 sin 2 9/ P. Sub-figures (l,2)-(2,3) show Sp ~ £ p . The equipotential surface, 
the colatitudinal components of the gravitational acceleration vector and the flux vector are 



-37- 



shown in the bottom panel, namely, sub- figures (3,1)- (3,3). 

Sub-figure (1,2) shows that pressure does not vary with colatitude 6 on the equipotential 
surface. It is the very feature that is required by the hydrostatic equilibrium on the surface. 
The numerical method of the solution to the 2D stellar structure equations presented in this 
paper is designed to achieve this feature. It is not trivial at all. 

Sub-figure (1,3) indicates that the presence of the magnetic flux loop beneath the surface 
affects the temperature distribution in site and above. This is reasonable since the thermal 
time scale near the base of the convection zone (where the loop is located) is much longer 
than the solar cycle so that the temperature perturbation travels little inwards in the cyclic 
period. In contrast, it can substantially travel outwards in short time since the thermal 
timescale above the torus field is very small. Another feature for the 2D temperature effect 
is that the temperature increases above the buried field. We see sunspots in the solar active 
regions. It is well-known that sunspots reduce the energy output of the Sun. We also know 
that the active regions increase the net energy output of the Sun as a whole. The idea 
that the temperature increase caused by the buried fields over-compensates the sunspot is a 
natural explanation to the net increase of the energy output in the active regions of the Sun. 

Sub-figures (2,1) and (3,1) are similar to each other. The distinction is their references: 
the former refers to the ID radius of the equipotential surface, and the latter refers to the 
equatorial radius. The maximal radius change takes place at the minimal (3 parameter. Both 
of them show the equipotential surface profile. 

Comparing sub-figure (2,3) with (1,1) we can see that the density change inversely 
follows the plasma (3 parameter and is of the same order of magnitude as 1/(3, which is 
in agreement with the analytical result: (p — < p >)/ < p >= 1/(1 + 1//3) ~ — 1//3. The 
sub- figure also shows that the density decrease maximizes in the loop. This will give rise to a 
buoyant force on the loop in the radial direction. Its component on the plane that is parallel 
to the equator plane cancels out since the loop is azimuthally symmetric. Its component in 
the meridional direction will generate an acceleration in the same direction, a m . Detailed 
calculation (see Appendix $B]) shows a m « 32 cm s~ 2 . The buoyant force is assumed to be 
balanced by the turbulent pressure generated by the down-flow plumes found in the realistic 
three-dimensional turbulent simulations of the solar convection zone near the surface of the 
Sun (e.g., Stein and Nordlund 1998; Robinson et al 2003). These simulations reveal that the 
up-flow and down-flow are not symmetric and the down-flow is stronger than the up-flow. 

In the real Sun, this condition is obeyed until the magnetic field reaches a critical value 
whereby the buoyancy forces dominate, magnetic loops making up the torus float up, produce 
magnetic activity in the solar surface, and the toroidal field is depleted. We do not model 
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these details in our code excepting in terms of the decrease of the toroidal field. 

The transverse components of the gravitational acceleration vector g and the flux F 
shown in sub-figures (3,2) and (3,3) are purely 2D effects. Their characteristics and other 
2D effects need to be investigated further and will be presented separately. 

6. CONCLUSIONS 

We present a new set of differential equations to describe the stellar equilibrium, in 
which two dimensional effects are explicitly taken into account. We improve the treatment 
presented in a previous paper of this series, by relaxing some approximations that had been 
made in that context; this task required one more differential equation, with the intro- 
duction of a new variable, i.e. the deviation of the radial component of gravity from the 
standard expression that is obtained when the Poisson equation is solved neglecting the 
angular derivatives. 

We have shown that by selecting an appropriate convergence criterion our code can 
reach the precision required by current and forthcoming observations. 

The code can now be used to test the effects of magnetic fields of any axisymmetric 
magnetic field configuration on the structure of the current Sun, and to investigate the 
change of the observable solar properties related to the variation of the magnetic field with 
the solar cycle. We have used the code to scan a very large region of the parameter space to 
test the code, and will present our findings in a separate paper. 

Finally, we wish to emphasize that because we are interested in modeling the effects of 
a dynamo-type field on the detailed envelope structure and global properties of the Sun, the 
code has been optimized for short timescales phenomena (down to 1 yr). Consequently, the 
time dependence of the code has so far been tested exclusively to address such problems, 
and we can not assume that the code could be used to model long term stellar evolution 
without further modifications. 

We want to acknowledge the following support for this work: LLH by NSF Grant ATM 
073770, and the Vetlesen Foundation; SS by the Vetlesen and the Brinson Foundations; 
SB by NSF grants ATM 0348837 and ATM 0737770; SLB by MSTC grant 2007CB8 15406, 
NSFC grants 10433030, 10773003, 10778601, and PD by NASA grant NAG5-13299. 
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A. COEFFICIENT MATRIX ELEMENTS 

Eq. (I98al) consists of a set of non-homogeneous linear algebraic equations. We work out 
these nonzero elements in this appendix. 



A.l. Useful Partial Derivatives 

The partial derivatives of the differential equations are required for the linearization. By 
defining the shorthand notation dxY = dY/dlogX, we can calculate the useful derivatives 
as follows. 

In fact, we need to calculate all the derivatives of V, T, W [i = 1, 2, 3, 4, 5), 71, and D 
(i = 1, 2) with respect to P', T', r', L, and U, respectively. For the sake of completeness 
and conciseness, we write down all nonzero partial derivatives and formulas except for the 
same as in Paper I. The derivatives of V, T, and C l are the same as in Paper I, where £} is 
equivalent to C in Paper I. 

The derivatives of 1Z may be nonzero only for k = j and 1 = i - 1, i. The unique nonzero 
derivative is 

d R K=-- K, 

which is different from Paper I. 

The derivatives of C e (£ = 2, 3) may be nonzero not only for k = j and 1 = i - 1, i. For 
the sake of simplicity, we rewrite Fg as follows: 

F = (F 1 + F 2 + F 3 )P, 

where 

F 2 
F 
P 



4acT 4 V 
3/tp 

1 pC P Tl m v conv V 

2 1+Wconv/fo 

1 pC pTl m v coav V &d 

2 1 + v conv /v 
GmQA 



Attt^Pt ' 

In order to obtain the nonzero derivatives of Fq, we also need the following formulas: 
OpF 1 = -F 1 (Kp + a-V P ) 
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r\ 7 — r 1 

o T F 


= —F (kt — S — Vt - 4) 


<JRr 


— V R 


d T F l 


= FVr 

x v Lj 


d P F 2 


= F 2 (a + Vp + Cpp) - - 

1 


d T F 2 


= F 2 (l -5 + Vt + Cpt) 


8 D F 2 


— -T V R 


(J LP 


F 2 T7 r 

— r V L 


d P F 3 


= F 3 (a + V P + C PP ) - - 


d T F 3 


= F 3 (l-5 + V T + C PT ) 


d R F 3 


= 


dpP 


= -P 


d T P 


= 


d R P 


= -AP 



'^I^F 2 (2a + C P p + K P ) 
Vco ™/ Vo F 2 (-25 + C PT + «r - 3) 

-+^convM) 



F 3 (2a + Cpp + /tp) 



1 +^conv/^0 

^conv/^0 



+ v conv /v 



F 3 (-25 + Cp T + kt-3) 



Here K p ee (Jg*^, Kr ee (fgf)^, C PP ee (§£g) r> C PT ee and *; = 

QacT 3 /p 2 C P l m K, V P = (91nV; d /91nP T ) T , and V r = (9 In V; d /<91nT)p T . As a result, we 



have 



3 

4 



where 



dpF e = P^dpF 1 + d P Pj2 F 

3 3 

dpFg = PJ2 d TF e + d T PJ2 F£ 

l=\ l=\ 

2 3 

d R F e = Pj2d R F e + dRpJ2 Fe 

e=i e=i 

2 

8 L F e = PJ2d L F e 
e=i 

J= = dpF e + d T F e -V + d R F 6 -^ + d L F e -C/V 

dr' A-KrlPp 
dP* ~ ~ GmQx ' 
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These finish the expressions for C 2 and £ 3 , and their derivatives: 

d P C 2 = C\F^d P F e -a) 

d T C 2 = C 2 (F~ l d T F e + 5) 

d R C 2 = C 2 {F e l d R F e -l) 

d L C 2 = C 2 F^d L F e 
hi 

hi 



d P C 3 = -a-C 3 



d T C 3 = 5-C 3 



d R C 3 = -C 3 



After all nonzero components and their derivatives are calculated, we can sum them to 
obtain 



3 



e=i 

3 

d P C = J2 dpjC - e 



e=i 

3 



d T c = J2 dTjC * 



3 

I 



OrC = drC 2 



hi 



A. 2. INTERIOR POINTS 

The interior points can be grouped into four blocks: 



Block I, 1 = i - 1 and k = j, 
Block II, 1 = i and k = j. 



A. 2.1. w = P 



For block I, 

dFp 

dFp 

dLi-ij 
dFp 

dPU 
dFp 

U1 i-ij 

For block II, 

dF'p 

dFp 
dLij 
dF P 
~dPJ 
dF P 
dT F - 



For block I, 

dF l 4 

dFi 

dLi-ij 

dPU 
dFif 1 

OI i-ij 



= 

= -^AsidpVi-i - 1 
= 

= 

= --AsidpVi + 1 
= 

A. 2.2. w = T 

= —AsidpTi^j 

= -^Asid L %_ X j 

= —AsidpTi^j 

= —AsidrTi^j - 1 
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For block II, 



For block I, 



For block II, 



dR 2 " J 



dLi_ij 

dF R 

0U t - 



--AsidLTij 



dF% _1 

"2 

-rr^r = — AsidpTa 

Qpi 2 J 

<9F ij 1 

— — y- = — AsidpTa + 1 

97?- 2 3 



A. 2. 3. w=R 



dF R 1 



dRU 2 



--AsidnKi-! - 1 



--Asj9p7vli_i 



1 



dTU 3 2 



" - --AsidnKi + 1 



^ 2' 

dL i:j 



= 



—4 = — AsidpTZi 
dP[ 2 

(9i^* 1 
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For block I, 



For block II, 



A. 2.4. w = L 



--AsidnCi-ij 



= 



i-lj 



dPU = 2 As ^ C ^ 

d F ij 1 

--Asidrd-u - 1 



dF l l _1 

~m ~ ~2 
= 



--Asid R C i: j 



= 



dL i:j 

dUij 

dF ij 1 

-^kr = — AsidpCa 

dP{ 2 3 

8F ij 1 

TT-r = — A Si 9 P A,- + l 

97V- 2 * J 



A.3. BOUNDARY POINTS 

A. 3.1. Center: w = R 



Central boundary points have only block II for w = R: 



dR[ 
dL ±j 



-45 - 



dFl 

— - = 

dU l3 

dFl 1 
dPj ~ 3 ao1 
dFk _ 1 
df[ 3 ~ 3 01 



A. 3. 2. Center: w = L 



Central boundary points have block II for w = L: 

OF} 3 



dFl 
dF^ 

dJl 

dP[ 
dF x L 3 



= 1 

= 

= -dpCij 

= -d T Cij 



A. 3. 3. Surface: w = R 
Surface boundary points have block I for w = R: 

1 



dF M+i 



dR' M 
dL M j 
dU Mj 

f)T' 

OI Mj 



= 

= 

= -Oi 

= — a 2 
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A. 3. 4- Surface: w 
Surface boundary points have block I for w = L: 

dF M+lj 

/ = 

9R Mj 

dF M+li 



dL 



Mj 



9F M+lj 

—± = 

dU M j 

dF M+lj 

-QpT- - ~L Mj a 4 



Mj 



-LiMjOj*, 



B. Buoyant acceleration of a magnetic flux loop in the meridional direction 

The magnetic flux loop used in this paper is assumed to be axisymmetric with respect to 
the polar axis. Its buoyant force (fg) is radial and can be decomposed into two components. 
One is parallel to the equatorial plane (/ e ), and the other is perpendicular to it (f m ). The 
former is canceled out since the loop is axisymmetric with respect to the polar axis (i.e., 
f e = 0), and the latter is in the meridional direction (f m ^ 0). In order to compute the 
buoyant acceleration of the loop in the meridional direction (a m = / m /mt), we have to 
compute f m and the mass of the loop mi- 

We must first calculate the boundary of the loop. The polar axis is assumed to be 
the z-axis. The equation for a torus azimuthally symmetric about the z-axis in Cartesian 
coordinates is 

( C _ v^f + ^-^^a 2 , (Bl) 

where c is the radius from the center of the hole to the center of the torus tube, a is the 
radius of the tube, and (0, 0, zq) is the center point coordinate of the hole. In the xz-plane 
the torus becomes two circles. One of them is 

(c - x) 2 + (z - z ) 2 = a 2 (B2) 

in Cartesian coordinates. We need to determine its boundary. In the spherical polar coordi- 
nates (r,9,4>), Eq. flB2j) becomes 

(c - r sin 9) 2 + {rcos6-c cot 6 ) 2 = a 2 , (B3) 
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where 6$ is the colatitude of the center of the circle. The radius range of the circle for each 
9 is given by the solutions for r of Eq. (|B3]): r_ < r < r + +, where r+ are defined by 

r± = c(sin6> + cos0cot0 o ) ± c[(sin6> + cos 9 cot 9 Q ) 2 - 1 - cot 2 9 + a 2 /c 2 ] 1/2 . (B4) 

The colatitude range of the circle for each radius r is determined by the solutions of Eq. (1B3[) 
for 9: 

b sin 29 ± [b 2 / sin 2 29 - 4(6 2 - 1) sin 2 fl ] 1/2 " 



#± = arccos 
where 



2 



(B5) 



= cV sin 2 9 + r 2 -a 2 ^ 

2cr v ; 



Since 9- >9 + , the boundary of Eq. (1B2I) can be expressed by 

C : r_ < r < r+, and 0+ < < 0_, and < < 2tt. (B7) 

We have two ways to define a torus field. One is to use the step function: Q = Qq within 
the loop confined by C, but Q = outside the loop, where Qo is a constant. The other way 
is to use the Gaussian profile to smooth the step function: Q = Q exp[— \ ^ e ~^ ], where 
a =1(9+ -9-). 

We can then express the meridional buoyant force component f m and mass in the loop 
in terms of the following integrals: 

f m = 2tcccos9 / rg(< p > —p)drd9, (B8) 
Jc 

rriL = 2ttc / rpdrd9. (B9) 



c 



The acceleration equals a B = f m /m L . Here < p > is the averaged density over the colatitude 
9 from to vr/2. 
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Convergence Criterion 



Fig. 1. — The numerical accuracy of the 2D solar models with a zero-rotation rate as a 
function of convergence criteria. The symbols mark the data points. 
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Fig. 2. — Maximal rotation effect as a function of the rotation rate Q. The symbols mark 
the data points, and the dashed line shows the relative error estimate. 
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Fig. 3. — Maximal rotation effect in all zones as a function of both the rotation rate Q and 
radius R/R & . 
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Fig. 4. — Contours of the 2D solar model with a uniform rotation rate Q 
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The top five sub-figures show the detail dependence of £p to £ p on R/R & and 0. The last 
fore sub-figures show the equipotential surface function x — 1, the transverse component of 
energy flux F, Fg/F Q , the radial perturbation component, and transverse component of the 
gravitational acceleration. 
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Fig. 5. — Maximal rotation effect as a function of the rotation rate Q. The symbols mark 
the data points, and the dashed line shows the relative error estimate. 
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Fig. 6. — Maximal rotation effect in all zones as a function of both the rotation rate Q and 
radius R/R & . 
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Fig. 7. — Contours of the 2D solar model with a uniform rotation-like toroidal magnetic 
field B = (0,0, (Anpy^Qr sin 9), where Q = 10" 6 s -1 . The top five sub-fi erures show the 
detail dependence of Sp to £ p on R/R Q and 9. The last fore sub-figures show the equipo- 
tential surface function x — 1, the transverse component of energy flux F, Fq/F q , the radial 
perturbation component, and transverse component of the gravitational acceleration. 
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Fig. 8. — Contours of the 2D solar variability model with a torus field, in which the applied 
magnetic field (measured in the plasma (3 parameter), the relative changes of the stellar 
structure variables (pressure, temperature, radius, luminosity and density) and the transverse 
components of the gravitational acceleration and flux vectors. 



